%% Load the data - Select the Files: 'Figure3andSupp'
close all; clear all; clc;
cd('.\Figure3andSupp')

%% Figure 3a and Extended data Figure 7
% Heat maps of theta in time (experiment)
bins = 10;
t = 6;
bin_edges = linspace(-pi, pi, bins + 1);
spacing = 5;
%---9 cp
load('./9cp/thetaAndtracks.mat')

[counts, centers] = hist(theta_norm',(bin_edges(2:end) + bin_edges(1:end-1))/2);
counts = (counts./sum(counts))./(2*pi/bins);
len = (t*30)/spacing;  %this is where you set the length in terms of frames

figure; imagesc([1 len*spacing]./30, centers, counts(:,1:len));
cmap = colormap('hot');
colormap(cmap); cbar = colorbar;
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({});
title('experiment');xlabel('t (s)'); ylabel('theta (radians)')
caxis([0 0.25])
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1.5;
cbar.LineWidth = 1;
cbar.Color = 'k';

print(['./9cp/theta_map_exp_t' num2str(t) '.pdf'],'-dpdf','-r1200')

% main figure
% writematrix("time (s)", 'Figure3.xlsx', 'Sheet', '3A-Top', 'Range', 'B1')
% writematrix(linspace(0, len*spacing./30, len), 'Figure3.xlsx', 'Sheet', '3A-Top', 'Range', 'C1')
% writematrix("orientation (rad)", 'Figure3.xlsx', 'Sheet', '3A-Top', 'Range', 'A2')
% writematrix(centers, 'Figure3.xlsx', 'Sheet', '3A-Top', 'Range', 'A3:A13')
% writematrix(counts(:,1:len), 'Figure3.xlsx', 'Sheet', '3A-Top', 'Range', 'C3')
% sup figure
% writematrix("time (s)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2D', 'Range', 'B1')
% writematrix(linspace(0, len*spacing./30, len), 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2D', 'Range', 'C1')
% writematrix("orientation (rad)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2D', 'Range', 'A2')
% writematrix(centers, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2D', 'Range', 'A3:A13')
% writematrix(counts(:,1:len), 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2D', 'Range', 'C3')


%---control

load('./control/thetaAndtracks.mat')

[counts, centers] = hist(theta_norm',(bin_edges(2:end) + bin_edges(1:end-1))/2);
counts = (counts./sum(counts))./(2*pi/bins);
len = (t*30)/spacing;  %this is where you set the length in terms of frames

figure; imagesc([1 len*spacing]./30, centers, counts(:,1:len)); colormap hot; colorbar
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({});
title('experiment');xlabel('t (s)'); ylabel('theta (radians)')
caxis([0 0.25])
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1.5;

print(['./control/theta_map_exp_t' num2str(t) '.pdf'],'-dpdf','-r1200')

% sup figure
% writematrix("time (s)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2A', 'Range', 'B1')
% writematrix(linspace(0, len*spacing./30, len), 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2A', 'Range', 'C1')
% writematrix("orientation (rad)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2A', 'Range', 'A2')
% writematrix(centers, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2A', 'Range', 'A3:A13')
% writematrix(counts(:,1:len), 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2A', 'Range', 'C3')
%----------------Heat maps of theta in time(simulation)-------------------
%
bins = 20;
t = 6;
%---9cp

load('./9cp/SimulationTracks.mat')
len = t/(.02*10); %each time step in simulation corresponds to 0.02 simulation time units (50 = 1 t)
[counts, centers] = hist(wrapToPi(Theta(1:len,:)+pi)',bins);
counts = counts./sum(counts);

temp = counts./counts(:,1); temp = (temp./sum(temp))./(2*pi/bins);
figure; imagesc(10*[0 len*.02], centers, temp);
caxis([0 0.25])

cmap = ColorMapDerek([0 0 0], [1 0 0], 'no');
cmap = colormap('hot');
colormap(cmap); cbar = colorbar;

title('simulation')
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({});
xlabel('$t = \tilde{t}W/u$ (s)','interpreter','latex'); ylabel('theta (radians)')
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth =1.5;

print(['./9cp/theta_map_sim_t' num2str(t) '.pdf'],'-dpdf','-r1200')

% main figure
% writematrix("time (s)", 'Figure3.xlsx', 'Sheet', '3A-Bottom', 'Range', 'B1')
% writematrix(linspace(0, 6, len), 'Figure3.xlsx', 'Sheet', '3A-Bottom', 'Range', 'C1')
% writematrix("orientation (rad)", 'Figure3.xlsx', 'Sheet', '3A-Bottom', 'Range', 'A2')
% writematrix(centers, 'Figure3.xlsx', 'Sheet', '3A-Bottom', 'Range', 'A3')
% writematrix(temp, 'Figure3.xlsx', 'Sheet', '3A-Bottom', 'Range', 'C3')
% sup figure
% writematrix("time (s)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2E', 'Range', 'B1')
% writematrix(linspace(0, 6, len), 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2E', 'Range', 'C1')
% writematrix("orientation (rad)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2E', 'Range', 'A2')
% writematrix(centers, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2E', 'Range', 'A3')
% writematrix(temp, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2E', 'Range', 'C3')
%
%---control

load('./control/SimulationTracks.mat')
len = t/(.02*10); %each time step in simulation corresponds to 0.02 simulation time units (50 = 1 t)
[counts, centers] = hist(wrapToPi(Theta(1:len,:)+pi)',bins);
counts = counts./sum(counts);

temp = counts./counts(:,1); temp = (temp./sum(temp))./(2*pi/bins);
figure; imagesc(10*[0 len*.02], centers, temp);
caxis([0 0.25])
colormap hot; colorbar; title('simulation')
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({});
xlabel('$t = \tilde{t}W/u$ (s)','interpreter','latex'); ylabel('theta (radians)')
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1.5;

print(['./control/theta_map_sim_t' num2str(t) '.pdf'],'-dpdf','-r1200')

% sup figure
% writematrix("time (s)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2B', 'Range', 'B1')
% writematrix(linspace(0, 6, len), 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2B', 'Range', 'C1')
% writematrix("orientation (rad)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2B', 'Range', 'A2')
% writematrix(centers, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2B', 'Range', 'A3')
% writematrix(temp, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2B', 'Range', 'C3')

%---9cp---theory
theta_x_theor = linspace(-pi, pi, 501);
t = 6;
tau = 14.52;
Drot = 0.069;
p = FPdistribution(theta_x_theor, 0:0.001:t, Drot, tau);

figure; imagesc(0:0.001:t, theta_x_theor, p');
caxis([0 0.25])
colormap hot; colorbar; title('theory')
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({});
xlabel('$t = \tilde{t}W/u$ (s)','interpreter','latex'); ylabel('theta (radians)')
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1.5;

print(['./9cp/theta_map_theory.pdf'],'-dpdf','-r1200')

% sup figure
% writematrix("time (s)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2F', 'Range', 'B1')
% writematrix([0:0.001:t], 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2F', 'Range', 'C1')
% writematrix("orientation (rad)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2F', 'Range', 'A2')
% writematrix(theta_x_theor', 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2F', 'Range', 'A3')
% writematrix(p', 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2F', 'Range', 'C3')

%---control---theory
theta_x_theor = linspace(-pi, pi, 501);
t = 6;
tau = 1000;
Drot = 0.069;
p = FPdistribution(theta_x_theor, 0:0.001:t, Drot, tau);

figure; imagesc(0:0.001:t, theta_x_theor, p');
caxis([0 0.25])
colormap hot; colorbar; title('theory')
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({});
xlabel('$t = \tilde{t}W/u$ (s)','interpreter','latex'); ylabel('theta (radians)')
ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1.5;

print(['./control/theta_map_theory.pdf'],'-dpdf','-r1200')

% sup figure
% writematrix("time (s)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2C', 'Range', 'B1')
% writematrix(0:0.001:t, 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2C', 'Range', 'C1')
% writematrix("orientation (rad)", 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2C', 'Range', 'A2')
% writematrix(theta_x_theor', 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2C', 'Range', 'A3')
% writematrix(p', 'ExtDataFigure2.xlsx', 'Sheet', 'Ext2C', 'Range', 'C3')
%% 3B
bins = 10;
bin_edges = linspace(-pi, pi, bins + 1);
sim_bins = 20;
t = 6; %real time
tau = 14.52; %rotation time of cell
Drot = 0.069;

spacing = 5;
theta_x_sim = linspace(-pi, pi, sim_bins);
theta_x_exp = linspace(-pi, pi, bins);
theta_x_theor = linspace(-pi, pi, 501);
theor_prof = (1/(2*pi))*(1+tan(theta_x_theor/2).^2)./(exp(-t/tau)+exp(t/tau)*tan(theta_x_theor/2).^2);

p = FPdistribution(theta_x_theor, 1:0.001:t, Drot, tau);
theor_prof_FP = p(end,:);
theor_prof_FP_steady = exp((1/(tau*Drot))*cos(theta_x_theor))/(2*pi*besseli(0,1/(tau*Drot)));
%---9cp---sim

load('./9cp/SimulationTracks.mat')
len = 50; %each time step in simulation corresponds to 0.02 simulation time units (50 = 1 t)
[counts, ~] = hist(wrapToPi(Theta(1:len,:)+pi)',sim_bins);
temp = counts./counts(:,1); temp = (temp./sum(temp))./(2*pi/(sim_bins));

sim_prof = temp(:,floor(t/(0.02*10))+1);

%---9cp---exp
load('./9cp/thetaAndtracks.mat')

[counts, ~] = hist(theta_norm',(bin_edges(2:end) + bin_edges(1:end-1))/2);
counts = (counts./(sum(counts)))./(2*pi/bins);

exp_prof = counts(:,floor(t*30/5));

figure;
plts(1) = plot(theta_x_theor, theor_prof_FP, 'k', 'linew', 2);hold on
plts(2) = plot(theta_x_sim, sim_prof,'r', 'linew', 2);
plts(3) = plot(theta_x_exp, exp_prof,'b', 'linew', 2);


xlim([-pi pi]);
ylim([0 0.25]);
xticks([-pi -pi/2 0 pi/2 pi])

ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1.5;
daspect([8 1 1])


lgd = legend([plts(3) plts(2) plts(1)], {'experiment','simulation','theory'}, 'fontsize', 16, 'location', 'south');
legend box off
print('./9cp/theta_cross_section.pdf','-dpdf','-r1200')

writematrix("theta exp (rad)", 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'A2')
writematrix(theta_x_exp, 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'B2')
writematrix("theta sim (rad)", 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'A5')
writematrix(theta_x_sim, 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'B5')
writematrix("theta theoretical (rad)", 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'A8')
writematrix(theta_x_theor, 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'B8')

writematrix("exp density (pdf)", 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'A3')
writematrix(exp_prof', 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'B3')
writematrix("sim density (pdf)", 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'A6')
writematrix(sim_prof', 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'B6')
writematrix("theoretical density (pdf)", 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'A9')
writematrix(theor_prof_FP, 'Figure3.xlsx', 'Sheet', '3B', 'Range', 'B9')

%% C

load('./SimData/Parameters.mat')
load('./SimData/MatchedValues.mat')
load('./SimData/RMS_diff(rho.rho_s).mat');
load('./ExpData/Turning_Rates(33-66).mat');

temp=fields(Parameters);
for i = 1:length(temp)
    eval([temp{i} '= Parameters.' temp{i} ';'])
end

%Exp densities
x = [49.56 139.65 229.73 319.82 409.91 500.00 590.09 680.18 770.27 860.35 950.44];
x = x./1000;
rho{1} = [0.001330 0.001127 0.001170 0.001086 0.001025 0.001044 0.001092 0.001109 0.001073 0.001131 0.001318];
rho{2} = [0.001496 0.001167 0.001036 0.001028 0.001023 0.00109 0.001085 0.001086 0.00108 0.001132 0.001355];
rho{3} = [0.001448 0.001166 0.001107 0.001046 0.001079 0.001076 0.00105 0.001036 0.001094 0.001117 0.001352];
rho{4} = [0.001404 0.001232 0.001182 0.001191 0.001149 0.001123 0.001025 0.000947 0.000959 0.001095 0.001275];
rho{5} = [0.000896 0.000834 0.000972 0.001045 0.001142 0.001102 0.001025 0.001064 0.001213 0.001477 0.001859];

%Exp Viscosity dist
Pos_Vis = [45.3 129.9 214.47 299.08 383.68 468.3 552.9 637.5 722.1 806.7 891.3]./1000;
Viscosity{1} = ones(1,length(Pos_Vis));
Viscosity{2} = [1.42 1.37 1.38 1.33 1.27 1.18 1.10 1.08 1.07 1.05 1.02];
Viscosity{3} = [2.26 2.17 2.09 1.96 1.81 1.60 1.41 1.24 1.16 1.11 1.10];
Viscosity{4} = [3.31 3.19 2.99 2.69 2.51 2.20 1.86 1.58 1.39 1.23 1.19];
Viscosity{5} = [5.57 5.47 5.29 4.75 4.17 3.58 3.01 2.28 1.75 1.45 1.26];

%Exp velocity profiles
% velocity dist %
V_exp{1} = [89.36 95.25 91.04 92.53 94.96 91.71 89.52 91.19 93.27 94.89 90.34];
V_exp{2} = [59.57 70.75 72.94 72.90 73.80 73.43 76.89 78.30 79.76 81.71 75.08];
V_exp{3} = [42.88 49.75 52.10 55.21 58.20 62.43 69.45 75.81 78.90 82.67 76.96];
V_exp{4} = [26.93 29.52 32.19 35.74 40.34 47.51 59.06 75.05 85.06 91.58 90.50];
V_exp{5} = [19.48 20.88 21.55 24.19 28.51 35.54 47.49 62.55 77.77 87.25 89.18];
x_vel = [49.56 139.65 229.73 319.82 409.91 500.00 590.09 680.18 770.27 860.35 950.44];
x_vel = x_vel./1000;

V_exp_std{1} = [13.08 10.90 13.62 11.73 10.53 12.14 12.86 13.56 11.20 9.98 10.23];
V_exp_std{2} = [0.84 1.50 2.24 1.74 1.99 1.55 0.93 1.48 1.41 0.36 1.79];
V_exp_std{3} = [3.66 4.85 4.84 5.31 5.85 4.46 3.26 2.27 3.60 1.59 3.73];
V_exp_std{4} = [1.47 2.12 1.98 2.01 2.05 1.63 3.13 5.65 7.83 9.9 9.88];
V_exp_std{5} = [0.90 1.06 1.04 1.43 2.22 2.08 2.92 3.76 4.00 3.60 3.26];

V_exp_temp = cellfun(@(z) flip(z), V_exp,'un',0);
V_exp_sigm={};
for i = 1:length(V_exp_temp)
    if i == 1
        V_exp_sigm{i} = @(z) 1+0.*z;
        slope_sigm(1) = 0;
        continue
    end
    param = sigm_fit(x_vel(2:end-1),V_exp_temp{i}(2:end-1),[],[],0);
    temp = @(x) (param(1)+(param(2)-param(1))./(1+10.^((param(3)-x)*param(4))));
    temp = @(x) temp(x)/temp(0);
    V_exp_sigm{i} = temp;
    temp = max(abs(diff(V_exp_sigm{i}(0:0.01:1))./0.01));
    slope_sigm(i) = temp;
end

% grad_eta = [0, 0.79, 2.2, 3.4, 7.2];
for i = 1:length(Viscosity)
    temp = interp1(Pos_Vis, Viscosity{i}, [0.33 0.66]);
    grad_eta_third(i)= (temp(1)-temp(2))./0.333;
    grad_eta_full(i) = Viscosity{i}(1) - Viscosity{i}(end);
    grad_eta(i) = temp(1)-temp(2);
end
scaling_arg = nanmean(grad_eta_third./grad_eta_full);

clim = 0.3;

x = [54.77 144.32 233.86 323.41 412.95 502.50 592.05 681.59 771.14 860.68 950.23];
x = x./1000;
accum_exp = cellfun(@(z) (2*trapz(x(1:6),z(1:6))-trapz(x,z))./trapz(x,z), rho, 'un', 0);
accum_exp = [accum_exp{:}];

figure('color',[1 1 1]);
imagesc([(list_V_r(1)-list_V_r(end))./list_V_r(1) (list_V_r(1)-list_V_r(1))./list_V_r(1)],log10([1/list_Tau_visc(1) 1/list_Tau_visc(end)]),PM); hold on

% save phase matrix
writematrix("grad_v/v_max", 'Figure3.xlsx', 'Sheet', '3C_map', 'Range', 'B1')
writematrix(linspace((list_V_r(1)-list_V_r(1))./list_V_r(1),(list_V_r(1)-list_V_r(end))./list_V_r(1), size(PM,2)), 'Figure3.xlsx', 'Sheet', '3C_map', 'Range', 'C1')
writematrix("omega (s^-1)", 'Figure3.xlsx', 'Sheet', '3C_map', 'Range', 'A2')
writematrix(1./list_Tau_visc', 'Figure3.xlsx', 'Sheet', '3C_map', 'Range', 'A3')
writematrix(flip(PM, 2), 'Figure3.xlsx', 'Sheet', '3C_map', 'Range', 'C3')

ind = linspace(0,1,256);
color_formap = ([100 ; 100 ; 200]./256).*0.8;
cmap = color_formap + ind.*([1; 1 ;1]-color_formap);
cmap = cmap';
cmap = flip(cmap);
colormap(cmap);
cbar = colorbar;
cbar.LineWidth = 1;
cbar.Color = 'k';

colors = {[0 0 0], [0.588 0.286 0.612], [0.929 0.694 0.125], [0.466 0.674 0.188], [0.85 0.325 0.0980]};

Tau_visc_exp = mean_omega;
Tau_visc_error = std_omega;
V_r_error_temp = cellfun(@(x) sqrt(x(end).^2 + x(1).^2), V_exp_std, 'un', 0); V_r_error_temp = [V_r_error_temp{:}];
Tau_visc_min = Tau_visc_exp - Tau_visc_error; Tau_visc_max = Tau_visc_exp + Tau_visc_error;
ind = find(Tau_visc_min < 0); Tau_visc_min(ind) = 0.00001;

exponent = 0.934; %swimming speed scaling with viscosity
eta_over_eta_0 = 1.01:0.01:20;
eta_over_eta_0_exp = cellfun(@(x) x(1)/x(end), Viscosity, 'un',0);
eta_over_eta_0_exp = [eta_over_eta_0_exp{:}];
W=1000;
alpha = 33;
y = scaling_arg*(2*alpha)/W*(1-(2./(eta_over_eta_0+1)));
y_exp = scaling_arg*(2*alpha)/W*(1-(2./(eta_over_eta_0_exp+1)));
x = 1-1./eta_over_eta_0.^(exponent);
x_exp = 1-1./eta_over_eta_0_exp.^(exponent);

plot(x,log10(y),'-k', 'linew',2);
for i =1:length(x_exp)
    plot(x_exp(i),log10(y_exp(i)),'dk','markerfacecolor',colors{i})
end
yticks([-3 -2 -1 -0.06]);
yticklabels({'10^{-3}','10^{-2}','10^{-1}', '10^{0}'})

% scaling curve and points
writematrix("grad_v/v_max", 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'A2')
writematrix(x, 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'B2')
writematrix("omega (s^-1)", 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'A3')
writematrix(y, 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'B3')
writematrix("theoretical grad_v/v_max", 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'A5')
writematrix(x_exp, 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'B5')
writematrix("theoretical omega (s^-1)", 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'A6')
writematrix(y_exp, 'Figure3.xlsx', 'Sheet', '3C_scalingcurve', 'Range', 'B6')


V_rr = [];
for i = 1:length(Tau_visc_min)
    V_rr(i) = (V_exp_sigm{i}(0)-V_exp_sigm{i}(1))./V_exp_sigm{i}(0);
    V_r_temp = (V_exp{i}(end)-V_exp{i}(1))./V_exp{i}(end);
    V_r_error(i) = sqrt((V_r_error_temp(i)./(V_exp{i}(end)-V_exp{i}(1))).^2 + (V_exp_std{i}(end)./V_exp{i}(end)).^2).*V_r_temp;
    x1 = ones(1,100)*V_rr(i);
    y1 = linspace(log10(Tau_visc_min(i)),log10(Tau_visc_max(i)),100);
    y2 = ones(1,100)*log10(Tau_visc_exp(i));
    x2 = linspace(V_rr(i)-V_r_error(i),V_rr(i)+V_r_error(i),100);
    line(x1, y1, 'color', colors{i}, 'linewidth', 2)
    line(x2, y2, 'color', colors{i}, 'linewidth', 2)
end

% Experimental points with error
writematrix("grad_v/v_max", 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'A2')
writematrix(V_rr, 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'B2')
writematrix("omega (s^-1)", 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'A3')
writematrix(Tau_visc_exp, 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'B3')
writematrix("grad_v/v_max error", 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'A5')
writematrix(V_r_error, 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'B5')
writematrix("omega error (s^-1)", 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'A6')
writematrix(Tau_visc_error, 'Figure3.xlsx', 'Sheet', '3C_expwitherror', 'Range', 'B6')

for i=1:5
    plot(V_rr(i),log10(Tau_visc_exp(i)),'o','markersize',6,'linew',2,'color',colors{i}, 'markerfacecolor', 'w');
    Vr_match_new =  (100 - Vr_match_exact(i))./100;% (100 - [100, 81.3793, 56.5517, 31.7241, 22.41])./100;
    Tau_visc_match = 1./Tau_match(i); % 1./[1000 398.11 31.62 31.62 10];
    plot(Vr_match_new,log10(Tau_visc_match),'sk','markersize',8,'linew',2,'color',colors{i},'markerfacecolor','w');
end

% Matched Points
writematrix("grad_v/v_max", 'Figure3.xlsx', 'Sheet', '3C_matchedsims', 'Range', 'A2')
writematrix((100 - Vr_match_exact)./100, 'Figure3.xlsx', 'Sheet', '3C_matchedsims', 'Range', 'B2')
writematrix("omega (s^-1)", 'Figure3.xlsx', 'Sheet', '3C_matchedsims', 'Range', 'A3')
writematrix(1./Tau_match, 'Figure3.xlsx', 'Sheet', '3C_matchedsims', 'Range', 'B3')

xlabel('$\Delta v/v_\mathrm{max}$','Interpreter','Latex','Fontsize',16);
ylabel('$1/\tau_\mathrm{visc}$','Interpreter','Latex','Fontsize',16);
set(gca,'Ydir','normal')

daspect([1 2.9 1])

ax = gca;
ax.XColor = 'k';
ax.YColor = 'k';
ax.LineWidth = 1;
title("rho|rho_s")
print('./PhaseMap/PhaseMap_RMS(rho.rho_s).pdf','-dpdf','-r1200')
%% FUNCTION

function p = FPdistribution(theta, t, Dr, tau)
Nt = length(t);
Ntheta = length(theta);
dt = abs(t(2)-t(1));
dtheta = abs(theta(2)-theta(1));
p0 = 1/2/pi*ones(size(theta));

p = zeros(Nt,Ntheta);

for ii = 1:Nt
    
    if ii == 1
        p(ii,:) = p0;
    else
        p(ii,2:end-1) =  p(ii-1,2:end-1) ...
            + dt*(   Dr * (p(ii-1,3:end) - 2*p(ii-1,2:end-1) + p(ii-1,1:end-2))/dtheta^2 ...
            + (1/tau) * (p(ii-1,3:end).*sin(theta(3:end)) - p(ii-1,1:end-2).*sin(theta(1:end-2)))/2/dtheta);
        
        p(ii,1) =  p(ii-1,1) ...
            + dt*(   Dr * (p(ii-1,2) - 2*p(ii-1,1) + p(ii-1,end))/dtheta^2 ...
            + (1/tau) * (p(ii-1,2).*sin(theta(2)) - p(ii-1,end).*sin(theta(end)))/2/dtheta);
        p(ii,end) =  p(ii-1,end) ...
            + dt*(   Dr * (p(ii-1,1) - 2*p(ii-1,end) + p(ii-1,end-1))/dtheta^2 ...
            + (1/tau) * (p(ii-1,1).*sin(theta(1)) - p(ii-1,end-1).*sin(theta(end-1)))/2/dtheta);
    end
    
    
end
end
%% FUNCTION
function [param,stat]=sigm_fit(x,y,fixed_params,initial_params,plot_flag)
% Optimization of parameters of the sigmoid function
%
% Syntax:
%       [param]=sigm_fit(x,y)
%
%       that is the same that
%       [param]=sigm_fit(x,y,[],[],[])     % no fixed_params, automatic initial_params
%
%       [param]=sigm_fit(x,y,fixed_params)        % automatic initial_params
%       [param]=sigm_fit(x,y,[],initial_params)   % use it when the estimation is poor
%       [param]=sigm_fit(x,y,fixed_params,initial_params,plot_flag)
%
% param = [min, max, x50, slope]
%
% if fixed_params=[NaN, NaN , NaN , NaN]        % or fixed_params=[]
% optimization of "min", "max", "x50" and "slope" (default)
%
% if fixed_params=[0, 1 , NaN , NaN]
% optimization of x50 and slope of a sigmoid of ranging from 0 to 1
%
%
% Additional information in the second output, STAT
% [param,stat]=sigm_fit(x,y,fixed_params,initial_params,plot_flag)
%
%
% Example:
% %% generate data vectors (x and y)
% fsigm = @(param,xval) param(1)+(param(2)-param(1))./(1+10.^((param(3)-xval)*param(4)))
% param=[0 1 5 1];  % "min", "max", "x50", "slope"
% x=0:0.1:10;
% y=fsigm(param,x) + 0.1*randn(size(x));
%
% %% standard parameter estimation
% [estimated_params]=sigm_fit(x,y)
%
% %% parameter estimation with forced 0.5 fixed min
% [estimated_params]=sigm_fit(x,y,[0.5 NaN NaN NaN])
%
% %% parameter estimation without plotting
% [estimated_params]=sigm_fit(x,y,[],[],0)
%
%
% Doubts, bugs: rpavao@gmail.com
% Downloaded from http://www.mathworks.com/matlabcentral/fileexchange/42641-sigmoid-logistic-curve-fit

% warning off

x=x(:);
y=y(:);

if nargin<=1 %fail
    fprintf('');
    help sigm_fit
    return
end

automatic_initial_params=[quantile(y,0.05) quantile(y,0.95) NaN 1];
if sum(y==quantile(y,0.5))==0
    temp=x(y==quantile(y(2:end),0.5));
else
    temp=x(y==quantile(y,0.5));
end
automatic_initial_params(3)=temp(1);

if nargin==2 %simplest valid input
    fixed_params=NaN(1,4);
    initial_params=automatic_initial_params;
    plot_flag=1;
end
if nargin==3
    initial_params=automatic_initial_params;
    plot_flag=1;
end
if nargin==4
    plot_flag=1;
end

if exist('fixed_params','var')
    if isempty(fixed_params)
        fixed_params=NaN(1,4);
    end
end
if exist('initial_params','var')
    if isempty(initial_params)
        initial_params=automatic_initial_params;
    end
end
if exist('plot_flag','var')
    if isempty(plot_flag)
        plot_flag=1;
    end
end

%p(1)=min; p(2)=max-min; p(3)=x50; p(4)=slope como em Y=Bottom + (Top-Bottom)/(1+10^((LogEC50-X)*HillSlope))
%f = @(p,x) p(1) + (p(2)-p(1)) ./ (1 + 10.^((p(3)-x)*p(4)));

f_str='f = @(param,xval)';
free_param_count=0;
bool_vec=NaN(1,4);
for i=1:4;
    if isnan(fixed_params(i))
        free_param_count=free_param_count+1;
        f_str=[f_str ' param(' num2str(free_param_count) ')'];
        bool_vec(i)=1;
    else
        f_str=[f_str ' ' num2str(fixed_params(i))];
        bool_vec(i)=0;
    end
    if i==1; f_str=[f_str ' + (']; end
    if i==2;
        if isnan(fixed_params(1))
            f_str=[f_str '-param(1) )./ (   1 + 10.^( ('];
        else
            f_str=[f_str '-' num2str(fixed_params(1)) ')./ (1 + 10.^(('];
        end
    end
    if i==3; f_str=[f_str ' - xval ) *']; end
    if i==4; f_str=[f_str ' )   );']; end
end

eval(f_str)

[BETA,RESID,J,COVB,MSE] = nlinfit(x,y,f,initial_params(bool_vec==1));
stat.param=BETA';

% confidence interval of the parameters
stat.paramCI = nlparci(BETA,RESID,'Jacobian',J);

% confidence interval of the estimation
[stat.ypred,delta] = nlpredci(f,x,BETA,RESID,'Covar',COVB);
stat.ypredlowerCI = stat.ypred - delta;
stat.ypredupperCI = stat.ypred + delta;

% plot(x,y,'ko') % observed data
% hold on
% plot(x,ypred,'k','LineWidth',2)
% plot(x,[lower,upper],'r--','LineWidth',1.5)

free_param_count=0;
for i=1:4;
    if isnan(fixed_params(i))
        free_param_count=free_param_count+1;
        param(i)=BETA(free_param_count);
    else
        param(i)=fixed_params(i);
    end
end

if plot_flag==1
    x_vector=min(x):(max(x)-min(x))/100:max(x);
    figure;
    plot(x,y,'k.',x_vector,f(param(isnan(fixed_params)),x_vector),'r-')
    xlim([min(x) max(x)])
end
end
%% FUNCTION
function cmap = ColorMapDerek(start_col,end_col,invert)
    cmap = [linspace(start_col(1),end_col(1),256)', linspace(start_col(2),end_col(2),256)', linspace(start_col(3),end_col(3),256)'];
    if strcmp(invert,'yes'); cmap = 1-cmap; end
end